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Previous  attempts  to  derive  the  homogeneous  closed  form  solution  to  the  problem 
of  the  dynamic  motion  of  a  zero  bending  rigidity  cylinder  in  a  viscous  stream 
have  expressed  the  solution  as  an  infinite  series  involving  Bessel  functions  of 
complex  argument  and  order,  which  are  often  impractical  to  evaluate  because  of 
their  complexity.  Moreover,  when  these  solutions  are  extended  to  nonhcrao^eneous 
situations,  a  harmonic  time  dependence  is  assumed  that  requires  ''’forcing'?  the 
system  by  an  arbitrary  time  function  using  multiple  solutions  combined  in  the 


do  , 1473 


20.  (Cont'd) 


r  Fourier  sense.  This  paper  presents  a  general  purpose  numerical  treatment 
formulated  to  overcome  these  difficulties.  The  numerical  approach  is  based 
on  finite  difference  schemes  applied  in  conjunction  with  powerful  numerical 
ordinary  differential  equation  methods.  The  theory  is  examined  with  respect 
to  consistency,  stability,  and  convergence  of  these  numerical  procedures. 

A  numerical  example  is  included  to  demonstrate  the  validity  of  the  treatment 
Although  an  explicit  boundary  condition  is  absent  from  this  study,  a  derived 
boundary  condition  is  demonstrated  to  be  adequate.^ 
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A  NUMERICAL  TREATMENT  OF  THE  DYNAMIC  MOTION  OF  A 
ZERO  BENDING  RIGIDITY  CYLINDER  IN  A  VISCOUS  STREAM 


1.  INTRODUCTION 

Paidoussls*  worked  out  a  solution  to  the  dynamic  motion  problem. 

2 

Ortloff  and  Ives  studied  a  special  case  of  the  same  problem  and  expressed 
their  solution  In  the  form  of  an  Infinite  series  Involving  Gamma  and  Bessel 
functions.  Both  the  orders  and  the  arguments  of  Bessel  functions  are 
generally  complex  and  can  be  large  In  magnitude.  Furthermore,  evaluation  of  a 
Bessel  function  of  complex  order  Is  laborious  and  time-consuming,  and  accuracy 
cannot  be  assured.  When  the  solution  proposed  by  Ortloff  and  Ives  is  applied 
to  the  nonhomogenous  problem  where  the  "upstream"  end  of  the  cylinder  is 
forced,  a  harmonic  time  dependence  is  assumed;  this  means  that  "forcing"  the 
system  by  an  arbitrary  time  function  will  require  multiple  solutions  combined 
In  the  Fourier  sense. 

To  overcome  these  difficulties,  a  general  purpose  numerical  approach  is 
Introduced.  This  approach  discretizes  ,  and  by  backward  and 

central  differences.  This  discretization  brings  the  dynamic  motion  equation 
Into  a  system  of  second  order  ordinary  differential  equations.  This  system  is 
decomposed  Into  a  system  of  first  order  ordinary  differential  equations.  A 
feasible  numerical  ordinary  differential  equation  method  is  then  used  to  solve 
this  system  with  optimal  efficiency. 

There  are  many  advantages  to  using  a  numerical  method  to  solve  the  problem 
of  dynamic  motion.  The  theory  Is  well  developed  with  respect  to  consistency, 
stabllvcy,  and  convergence.  Numerical  methods  are  systematic  to  Implement, 


and  effective  techniques  can  be  used  to  accurately  accelerate  computations. 
When  a  numerical  approach  Is  used,  the  laborious  evaluation  of  special 
functions  Is  bypassed,  maximizing  accuracy  and  efficiency. 

This  report  begins  with  a  description  of  the  dynamic  motion  problem  and 
the  associated  Initial  and  boundary  conditions.  A  numerical  approach  Is 
Introduced  and  the  supporting  theory  and  mathematical  formulation  are 
discussed.  An  example  Is  given  to  demonstrate  the  validity  of  our  numerical 
solution  to  a  well  posed  dynamic  motion  problem.  The  computer  programs  are 
Included  In  the  appendix. 
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2.  POSED  PROBLEMS 


The  motion  of  a  wire  suspended  in  a  fluid  stream,  considered  by  Ortloff 

p 

and  Ives,  can  be  described  mathematically  by  the  partial  differential 
equation. 


El  -2-jf  +  (M  +  m)  +  MU2  +  2MU 

3X**  3 1  3X^  3t3X 


3 

3X 


'T  M  i/2  n  v)  M. 

r  is  u  (Ux)  it 


+  7cnSu  3t  +  U« 


0, 


(2.1) 


where 

El  »  bending  rigidity, 

M  >  lateral  virtual  mass  of  fluid  per  unit  length  of  wire  accelerated  by 
the  accelerating  wire, 
m  *  mass  of  the  wire  per  unit  length, 

U  »  velocity  of  the  free  stream, 

Cj  «  drag  coefficient  due  to  pressure  acting  on  the  wire  surface, 

D  «  wire  diameter, 

L  a  total  wire  length,  and 

CN  a  drag  coefficient  due  to  shear  forces  acting  on  the  wire  surface. 
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The  special  case  of  a  zero  bending  rigidity  (or  an  infinitely  flexible) 
cylinder  is  realized  by  setting  El  =  0.  To  express  the  above  equation  in 
dimensionless  terms,  use 

r  =  ( t/L )U  ,  e  =  M/(M+m)  ,  i  =  x/L  , 

c  =  L/D  ,  n  =  y/L  . 

The  above  equation  becomes 


,2  ,2  r  ,  i  .  (CT  +  CM) 

77  77  L  7  T  c  (1  "e)J  Tt  - 2 -  eB 


CH«  I?-0- 


(2.2) 


The  associated  initial  boundary  conditions  are  described  by 


n  ■  0 

4» 

0 

(fixed  end  condition); 

(2.3)a 

|n|  is  finite. 

4* 

1 

(bounded  free  end  deflection); 

(2.3)b 

n  -  nj  (4) 

r  a 

0 

(prescribed  initial  deflection); 

(2.3)c 

and 


Jt-q  =0,  0  (zero  initial  velocity).  ( 2 . 3 ) d 

Ortloff  and  Ives  solved  the  problem  posed  by  equation  (2.2)  using 
conditions  described  in  equation  (2.3).  Their  solution  is  expressed  in  terms 
of  Bessel  functions. 
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The  initial  boundary  value  problem,  equation  (2.3),  for  the  partial 
differential  equation  (2.2)  is  said  to  be  well  posed  in  the  sense  of 
HadamardJ  if  and  only  if  its  solution  exists,  is  unique,  and  depends 
continuously  on  the  data  assigned.  After  the  problem  is  formulated  using 
finite  difference  and  ordinary  differential  equations,  it  will  be  seen  that 
the  problem  is  well  posed.  We  will  seek  a  unique  solution  by  means  of  the 
numerical  techniques  presented  in  the  next  section.  When  the  boundary 
conditions  become  uncertain,  there  Is  not  enough  Information  available  to 
solve  equation  (2.2);  we  term  this  problem  ill  posed.  However,  a  derived 
boundary  condition  is  developed,  which  is  shown  to  be  adequate  for  our  problem. 

3.  THE  NUMERICAL  TREATMENT 

In  search  for  a  general  purpose,  accurate  solution  to  the  well  posed 
problem  (2.2),  subject  to  conditions  described  in  equation  (2.3),  the  method 
of  attack  is  to  discretize  ,  n^T  ,  and  by  central  and  backward  finite 
differences  and  then  to  transform  equation  (2.2)  into  a  system  of  second  order 
ordinary  differential  equations  (known  as  the  method  of  lines4).  We 
discovered  that  Generalized  Adams  Bashforth  (GAB)  methods^  can  be  used  to 
solve  this  system  efficiently. 

Expressing  equation  (2.2)  in  short  form  and  writing  u  as  n  gives 


urr  +  a(E)  u^  ♦  bu^  +  2s(u?)r  ♦  cur  -  0,  (3.1) 
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where 

am  -  a  [i  4^(1-*)]  * 

b  “2^CT  +  CN^c8*  and 
c 


3.1  FINITE  DIFFERENCE  DISCRETIZATION 

Applying  the  second  order  central  and  backward  finite  difference 
discretization  to  equation  (3.1)  in  the  <  direction,  we  obtain 


(u  )  +  a($  )  — - — ^-  +  £(u  -  u  ,) 

'  m’rr  '*m;  ^2  h  v  m  m-17 


+  4^  (u  -  u  ,)  +  c(u  )  *  0, 

h  '  m  nv-l'r  '  m'r  * 


(3.2) 


where  h  «  for  index  m  «  1,  2,  ... 

A  simplification  of  equation  (3.2)  gives 


<um>„  *  Or  *  c)(u*)r  -  r  ( Vl>,  *  ^ u, 


(3.3) 
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and 


“m  ‘  Vl 


<“m>{  - 


(3.6) 


Substituting  the  power  expansions  of  u^  and  u^  into  equation  (3.2) 
and  using  equations  (3.S)  and  (3.6),  we  find  that  equation  (3.2)  is  then 
expressed  in  a  difference  operator  form. 


t  h[u;h]  «  +  a(t) 

\3f 


a2  .  h2  ,4 

7?  *7? 


♦b 


a  ha2  L  a  T  a  h2  a2 

5T-7^J*2‘jr[»?-r^ 


r[u]  -  rh  [u;h] 


k2  h  ,2  h2  a2 

[-«(*)  T7i7*7bi7t  28  5-^7)“' 
1  l£  ar  £  a r  ar' 


It  is  seen  that 


(3.7) 


lim 

h*o 


UCu]  -  £h 


[u;h])  *  0. 


Therefore,  the  difference  operator  is  consistent  with  the  true  operator  in  the 
sense  of  Keller7.  Thus,  the  consistency  requirement  is  established. 
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Now  that  we  have  a  consistent  difference  operator,  we  seek  stable 
numerical  ordinary  differential  equation  methods  to  solve  equation  (3.2). 

3.2  ORDINARY  DIFFERENTIAL  EQUATION  SOLUTION 

To  seek  the  solution  to  equation  (3.2),  refer  to  equation  (3.3). 

Write  au 

m  *  w_ 


Equation  (3.8)  is  a  set  of  equations  that  represent  equation  (3.3)  as  a 
system  of  first  order  ordinary  differential  equations.  For  illustration, 
using  m  *  1,  2,  we  can  obtain 
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which  is  in  the  form 

u'  -  A($)u  +  g($,  r,  u).  (3.10) 

The  elements  contained  in  the  components  of  the  g-vector  have  the  following 
meanings. 
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Uq  -  top  boundary  fixed  end  condition:  Uq  -  0,  4*0; 

Wq  «  Initial  condition:  uT  ■  0,  r*  0;  and 

-  bottom  boundary  bounded  free-end  deflextlon  condition: 

|u^|  is  finite  at  1. 

The  matrix  elements*  A^j,  of  matrix  A  can  be  determined  by  the  following 
setups  In  which  we  define  A^  -  0  If  1-j  <  0  for  j  -  1,  2,  3. 

When  Index  1  Is  odd*  A.  ^  .  1.  When  Index  1  Is  even* 


A1 . 1—3  *  h  » 

A  »  ll 
1,1-2  h  * 


where  a(6)  Is  evaluated  at  a(6J,  l»  1 

I  j. 

Now,  the  problem  Is  to  select  an  effective  numerical  ordinary  differential 
equation  method  to  solve  equation  (3.8).  A  close  examination  suggests  that 
the  Generalized  Adams-Bashforth  (GAB)  method  offers  an  efficient  solution.  In 


II 
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the  present  application,  because  a  low  order  GAB  method  can  do  the  job,  high 
order  GAB  methods  are  not  necessary;  hence,  first  order  GAB  methods  were 
developed  Into  computer  programs  In  FORTRAN  language.  However,  the  program 
package  Is  flexible  so  that  high  order  GAB  methods  can  be  incorporated  when 
required. 

We  Introduce  the  first  order  GAB 


n+1  Ah  n  ,  .  »  .  . 

u  -  e  u  +  h#^0(Ah)  gn 

to  solve  equation  (3.10),  where 

*l,0(Ah)  *  -W1  (I  -  «Ah). 


(3.11) 


(3.12) 


The  theory  with  respect  to  consistency,  stability,  and  convergence  has 
been  very  well  developed  for  Nonlinear  Multistep  (NLMS)  methods.8  The  GAB 
method  is  a  member  of  the  NLMS  family.  We  summarize  the  theory  below. 

NLMS  methods  take  the  expression. 


io  ®jeAh*k  ^  Vi  *  h  i£o  *ki  Vr 

3.2.1  Stability 

The  characteristic  polynomial  of  NLMS  Is  defined  by 


(3.13) 


-**Ah  iJoV1- 


(3.14) 


12 


Using  the  GAB  method,  the  selection  of  is  such  that 

°k  *  1*  °k-l  *  _1*  We  see  that  the  root  of  p(x,*)  has  unity  and  is 

simple;  therefore,  method  (3.11)  Is  stable. 

3.2.2  Consistency 

The  GAB  method,  equation  (3.14)  satisfies  the  consistency 

condition 

llo  liio  aie  Ah(k-i)Vi  -  hU0  #ki(Ah,Vl|  "  0  (3 

for  k  -  1,  «  1  and  -  -1.  Therefore,  GAB  method  is  consistent. 

3.2.3  Convergence 

According  to  the  convergence  theorem  of  NLMS  methods,  "A  stable  and 
consistent  NLMS  method  Is  convergent."  Therefore  the  GAB  method  applied  to 
problem  (3.10)  Is  a  convergent  method. 

4.  BOUNDARY  CONDITIONS 


In  real  applications,  at  1,  the  bounded  free  end  deflexion  boundary 
condition  is  expected  to  be  such  that  n(l,r)  is  finite.  However,  the 
appropriate  function  n(l,r)  to  be  used  for  the  boundary  condition  appears 
uncertain  in  reality.  This  lack  of  Information  defines  problem  (2.2)  as  an 
111  posed  problem.  For  general  partial  differential  equations,  it  is  always 
difficult  to  formulate  correct  conditions  leading  to  a  well  posed  problem. 
Problems  may  look  reasonable,  yet  cannot  be  solved.  It  is  hoped  that  the 
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bounded  free  end  deflexion  boundary  condition  may  be  obtained  through 
experimentation,  but  the  exact  mathematical  expression  for  n(l,0  must  still 
be  worked  out.  We  will  attempt  to  change  the  111  posed  problem  to  a  well 
posed  one  so  that  a  solution  exists  and  can  be  solved  by  the  numerical 
techniques  we  have  developed. 

In  the  theory  of  second  order  partial  differential  equations  there  exists 
a  class  of  well  posed  problems,  such  as  the  Cauchy  problem  for  wave  equations, 
the  Dlrlchlet  condition  for  Laplace  equations,  and  the  mixed  Initial  boundary 
value  problem  for  heat  equations.  Our  first  step  Is  to  examine  the  most 
general  boundary  conditions.  Let  uN  denote  the  normal  derivative.  The 
first  boundary  value  problem  of  the  Dlrlchlet  type  Indicates 

u  -  f  (4.1) 

on  the  boundary.  The  second  boundary  value  problem  of  the  Neumann  type 
Indicates 


uN  -  f  (4.2) 

on  the  boundary.  The  third  boundary  value  problem  of  the  mixed  type  Indicates 

Ujy  ♦  au  -  f  (4.3) 
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on  the  boundary.  Note  that  the  third  boundary  value  problem  is  well  posed 
only  for  the  restricted  choice  of  a.  We  will  assume  that  the  free  end 
deflexion  boundary  condition  takes  the  expression 

xuN  +  au  m  f.  (4.4) 


When 


x 

x 

x 


0,  a  »  1,  (4.4)  reduces  to  (4.1); 

1,  a  -  0,  (4.4)  reduces  to  (4.2);  and 
1,  a  arbitrary,  (4.4)  reduces  to  (4.3). 


In  our  application,  as  given  by  the  numerical  example  in  the  next  section, 
n(£,0N  -  n (4,t)  ,  x  -  0,  o  -  1  gives 

n(i,r)  -  f($,  )  and  Jf (1,  r  )|  is  finite.  This  gives  Ortloff's  and  Ives' 

bounded  free  end-  deflexion  boundary  condition. 

The  procedure  to  be  followed  here  for  determining  a  free  end  boundary 

condition  is  to  derive  an  approximate  boundary  condition  and  then  to  use  that 

boundary  condition  to  compare  the  solution  with  a  direct  application  of  the 

2 

Ortloff  and  Ives  solution.  We  develop  a  form  of  the  boundary  condition  for 
the  second  order  partial  differential  equation  by  following  the  approach  used 


15 


TR  6343 


by  Paldoussls  for  his  fourth  order  partial  differential  equation;  that  Is,  by 
Integrating  the  transverse  momentum  equation  over  a  short  tapered  end  which  Is 
attached  to  the  free  end  In  order  to  generate  the  required  boundary 
condition.  Paldoussls  assumed  that  the  cross  sectional  area  tapers  smoothly 
from  S  to  zero  In  a  distance  (4)  sufficiently  short  that  the  forces  acting  on 
the  tapered  end  can  be  lumped  and  considered  In  appropriate  boundary 
conditions,  for  our  present  problem  the  boundary  condition  Is  obtained  from 

'/  (h  *  <%r)  ["<«>»]  "*  ♦/  V*  -  /  fe  (T<«>&)  d* 

•'I-!  '  '  L-l  L-l 

L 

.2 

m(x)  dx  »  0 
at^ 

L-I  (4.5) 


where  the  parts  of  the  equation  express  rate  of  change  of  fluid  momentum, 
hydrodynamic  forces,  and  cylinder  Inertia,  respectively,  and  where  f  Is  a 
factor  Introduced  by  Paldoussls  to  account  for  the  Intractable  flow  conditions 
at  the  free  end  and  v  is  the  transverse  velocity  of  the  fluid  relative  to  the 
cylinder.  Therefore , 


(4.6) 
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FN  "  7  (5)  UCNV’ 

T(x)  -  T(L)  ♦£({})  U2Ct  (L-x); 

and  T(L)  is  a  consequence  of  form  drag  at  the  free  end. 

An  important  assumption  necessary  to  perform  the  integration  is  that  the 
length  of  the  tapered  section  (1)  is  small  enough  that  the  lateral  velocity 
(V)  may  be  considered  constant  over  1.  We  find 

fM!  ($*u4r)  (»  tu  ») 

*1  (S)  “V© 

-  T(L)  I  +  2  I  ix  ♦  0  (I2)  -  0 

ax^  d  ar  (4.9) 

for  x  ■  L,  all  t. 

After  nondimension  of  this  equation  as  before  and  neglecting  terms  of 
order  (1 2)  and  |  ,  we  have 

for  •  1,  all 


an 

a? 


(4.10) 


(4.7) 

(4.8) 
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On  physical  grounds  it  is  reasonable  to  neglect  Cy  relative  to 
g 

Cn/4*  making  the  final  boundary  condition 


3n  a.  an 

3(  3T 


0  for 


-  1,  all  t. 


(4.11) 


which  amounts  to  a  "radiation  condition";  that  is,  no  reflected  energy 
exists.  In  the  following  sections,  we  refer  to  boundary  condition  (4.11)  as 
Kennedy  boundary  condition. 


5.  A  NUMERICAL  EXAMPLE 


The  test  example  is  obtained  through  linearization  of  a  fourth  order 
nonlinear  cable  equation,10  which  is  given  by 


(5.1) 


where  m,  M,  U  take  the  same  definitions  as  given  in  section  2.  T,  CN  are 
defined  as 


T 


U  -  X), 


(5.2) 


OU2  L, 


(5.3) 
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1  1  1_  P 

CT  -  2  7  q2  1,  and 


"  7  S  CNU* 


(5.4) 

(5.5) 


where  CN  and  Cy  satisfy  the  definitions  given  in  section  2. 


Assuming  and  a_  commute,  expanding  equation  (5.1)  gives 
ax  at 


at*  (  at* 


♦  2U  +  U2 
atax 


17)  1?  ,x,x  M\st  *V  (5.6) 


Performing  aT  and  using  definitions  (5.2)  through  (5.6),  we  get 
ax 

^  nr +  ^  p-iSr  (L'X)]  +  «  w  {ct  +  cn>  +  2U  h/k  *  \  cn  S  It  *  °- 

(5.7) 


Equation  (5.7)  is  the  same  as  Ortloff’s  equation  (2.1)  before 
nondimensionalization.^ 

Select 


M  .  m  -  0.00273 
U  -  15  ft/sec 
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L  -  2000  ft 

CT  -  1.8 

CN  -  1.1259. 

Then, 

e  -  48,000 
B  -  1 

7. 

The  solution  to  equation  (5.7)  is  expressed  by 

y(t,  x)  -  eiwt  J7(x), 

where  7  is  approximately  21.89  and  J7(x)  initial  values  are  calculated  using  a 
UNIVAC  1108  Bessel  function  subroutine.^ 

The  fixed  end  boundary  condition  initially  is  zero.  The  free  end  boundary 
condition  uses  n(l,T)  »  e1WTJ7(l). 

This  problem  was  tested  again  using  Kennedy's  free  end  boundary 
condition.  Results  are  surprisingly  in  agreement  with  the  known  solution. 

The  test  results  seem  to  show  that  the  Kennedy  free  end  boundary  condition  is 
adequate.  Results  are  presented  in  graphic  form.  Two  sets  of  graphs  are 
given:  one  displays  M4»t)|  versus  t,  the  other  displays  the  real  {n(€,r)[ 
versus  ▼.  Both  plots  are  constructed  at  g  ■  0.2,  0.8. 
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S  »  0.2  {  =  0.2 


Figure  1:  Solution  magnitude  vs  time  Figure  2:  Real  part  solution  vs  time 


{  =  08  £  =  0.8 


Figure  3:  Solution  magnitude  vs  time  Figure  4:  Real  part  solution  vs  time 
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6.  CONCLUSIONS 

A  numerical  solution  to  the  dynamic  motion  of  a  zero  bending  rigidity 
cylinder  in  a  viscous  stream  has  been  introduced.  The  numerical  procedures 
developed  to  obtain  the  solution  are  theoretically  convergent  and 
computationally  accurate. 

For  given  appropriate  boundary  conditions  and  accurate  initial  values, 
this  model  will  produce  an  accurate  unique  solution.  For  uncertain  boundary 
conditions,  this  model  can  be  used  as  a  tool  to  study  the  boundary  effects  and 
possible  to  construct  the  ad  hoc  boundary  conditions. 
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APPENDIX 

COMPUTER  PROGRAMS  STRUCTURE  AND  COMPUTER  LISTING 

COMPUTER  PROGRAM  STRUCTURE 


MAIN 
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ACRONYMS 

MAIN  main  program  which  controls  the  setup  of  inputs  and  the  preparation  of 
outputs 

START  supplies  the  initial  values 

OIFEQ  controls  the  present  r-step  and  calls  for  NLMS(GAB)  method 

NLMS  1st  order  Generalized  Adams-Bashforth  method 

INVER  calls  for  complex  matrix  inversion 

CGJR  complex  matrix  inversion  using  Gauss-Jorden  reduction 

GFN  calculates  the  g-vector 

BC  fixed  end  and  free  end  boundary  conditions 

PACE  a  rational  function  approximation  for  matrix  exponentials 

INVERT  calls  for  double  precision  matrix  inversion 

DGJR  double  precision  matrix  inversion  using  Gauss-Jorden  reduction 

The  user  needs  to  deal  with  MAIN,  START,  GFN,  and  BC.  The  user  need  not 
be  concerned  with  DIFEQ,  NLMS,  INVER,  INVERT,  CGJR,  and  DGJR. 
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COMPUTER  PROGRAM  LISTING 

MAIN 


COMMON  PEFS.PCT.PCN.PBETA.FQ.PA.PB.PC.PD.DZ, IBND 
COMMON  AO.8, 18>,T<3> 

DIMENSION  Y < 2 *  1 8 ) , YZERO< IS) , YNEW< 18) » EXACT (IS) 

COMPLEX  A  ,  Y , YZERO , YNEW ,  XX  ,  OMGA , SAVE  .EXACT 
*****  THIS  PACKAGE  SOLVED  A  2ND  ORDER  P.D.E.  BY  THE  METHOD  OF  LINES  AN 

GENERALIZED  APAMS-BASHFORTH  METHODS 
*****  REFERENCE:  ORTLOFF  AND  IVES 

*****  input  PARAMETERS  HAVE  THE  FOLLOWING  DEFINITIONS: 

*****  N  =  NUMBER  OF  2ND  ORDER  ODE 

*****  TMAX  =  MAXIMUM  TAO 

*****  TINT  =  EVERY  TAO  INTERVAL  TO  BE  PRINTED  OUT 
*****  PXI  ■  AT  THIS  XI.  THE  OUTPUT  IS  REQUESTED 

*****  FQ  *  FREQUENCY.  NONDIMEN SIGNAL  OMEGA 

*****  IBND  *  BOUNDARY  CONDITION  INDICATOR 

*  1  BUILT-IN  KENNEDY  BOUNDARY  CONDITION 

*  2  USER-SUPPLIED  BOUNDARY  CONDITION 

*****  PEPS  *  EPSILON 
*****  PCT  *  C  SUB  T 
*****  PCN  *  C  SUB  N 
*****  PBETA  =  BETA 
*****  H  *  TAO  STEP  SIZE 

*****  READ  INPUTS  HERE 

READ<5,*)  N, TMAX. TINT. PXI. FQ» IBND 
RE AD (5.*)  PEPS. PCT, PCN, PBETA.H 
N— 2*<N— 1 ) 

DZ=1 »0/<FL0AT (N/2+1 ) ) 

PB-0 . 5*PBETA*PCT*PEPS 
PA»PBETA*<1.0-PB) 

PD-0 . 5*PCN*PEPS*PBETA 
PC-PB+PD 
C  *****  TO  SET-UP  MATRIX  A 
DO  20  I-l.N 
DO  20  J-l.N 

20  A<I,J)»CMPLX<0. 0,0.0) 

BH-PC/DZ 
TBH»2.*PBETA/DZ 
DG— (TBH+PB) 

DO  28  I-l.N 
IF( I  ,GE.  4)  GO  TO  26 
IF( I  . EQ.  2)  GO  TO  25 
A< I »  T  +  l )«CMPLX (1.0, 0.0) 

GO  TO  20 

25  J-I/2 

X-PBETA*  < 1 . - . 5*PCT*PEPS* ( 1 . - J*DZ  > ) 

A< I , 1-1 )"CMPLX<2»*X/<DZ*DZ)-BH,0»0) 

A  < I , I > -CMPLX  <  DG , 0 • 0  > 

A  <  I ,  I  + 1 ) -CMPLX  < -X/ <  DZ*DZ ) , 0 . 0 ) 

GO  TO  28 

26  IF<MOD< 1.2)  .EG.  0)  GO  TO  27 
A< 1 , 1+1 ) -CMPLX  < 1 .0,0.0) 

GO  TO  28 


27  J=I/2 

X=PB ETA*  ( 1 .  -  .  s*pct*pfps*  ( 1 .  -  J*BZ )  ) 

A<I,I~3>=CMPLX<BH,0*0> 

4 < I , 1-2 > =CMPLX < TBH , 0 . 0 > 

A  < I » 1-1 ) =CMPLX <  2 . *X/ <  BZ*DZ  > -BH ,0.0) 

A  ( I  » I  >  =CMPLX  <  DR ,0.0) 

IP (I  «EQ»  N >  GO  TO  28 

All *141  >=CMPLX<-X/<DZ*DZ>fO.O> 

28  CONTINUE 
T<1>*0.0 

C  *****  TO  0BT4IN  INITIAL  VALUES  FROM  'START' 

CALL  START (N.T.VZERO) 

TX-0 . 0 

SAVE=YZERO  <  N-l ) 

10  CONTINUE 
tx=tx+tint 

CALL  D I FEO <  H  *  N , TX , Y , YZERO , Z » YNEW , SAVE ) 

C  *****  RESULTS  IN  YNEW(I) 

C  *****  Z  CONTAINS  PRESENT  TAO 

C  *****  USER  USES  ABOVE  INFORMATION  FOR  HIS  PLOT 
SAVE=YNEW<N-1 > 

IPT*IPT+1 

IF< IPT  .LE.  4)  GO  TO  8 
IF <MOD< IPT » 1 )  .EG.  0)  GO  TO  8 
GO  TO  100 
8  CONTINUE 

WRITE<5» 1 )  HfZ 

1  FORMAT ( /10X  »  '  H  «' ,E15.8»5X, 'T  ='»E15.8/> 

2  FORMAT ( 3X,F8 .2,10X,E15.8» 10X »E15 • 8) 

C  *****  PRINT  OUT  EXACT  SOLUTION 

C  *****  THIS  PORTION  IS  FOR  TEST  EXAMPLE  ONLY 
CALL  START <N»T • EXACT) 

BO  100  1=1, N. 2 
NN=<I+1 )/2 
DLZ=NN*DZ 

WRITE(5,2)  DLZ,  YNEW<I) 

WRITE < 5,3)  EXACT < I ) 

3  FORMAT <21X,E15.8,10X,E15.8/) 

100  CONTINUE 

IF<TX  .GE.  TMAX)  STOP 
GO  TO  10 
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NLMS 

SUBROUTINE  NLMS ( H , v , N , YN , IS , SAUE ) 

COMMON  PEPS  , PCT , PCN  ,  PBETA , EG , PA  , PB * PC  , PD ,  DZ ,  I BND 
PARAMETER  KM=l8 
COMMON  A(KM,KM) ,T(3) 

DIMENSION  AH ( KM , KM ) , EAH ( KM , KM ) *  G ( KM > 

DIMENSION  PI (KM, KM) ,UNIT <KM,KM) 

DIMENSION  Y < 2 , KM ) , YN ( KM ) 

DIMENSION  FE(KM,KM) ,A1 <KM,KM) 

COMPLEX  A, AH, EAH, G, PI , PH, UNIT, Y, YN, SAVE 
COMPLEX  FE,A1 
DATA  IND/O/ 

IF< IND  ,GT «  0)  GO  TO  14 
DO  2  1  =  1, N 

DO  2  J=1,N 

2  AH ( I , J ) =H*A  < I , J ) 

IF(N-l)  7,8,7 

8  Al(l,l)=l./AH<1,1) 

GO  TO  9 

7  CALL  INUER(AH,N,A1 ) 

9  IF(N-I )  10,11,10 
11  EAH  (1,1)  =CEXP  <  AH  <  1. ,  1 )  ) 

GO  TO  14 

10  CALL  PADE  < A , H » E AH , N ) 

IND  =  IND  +  1 

14  IF (IS  *0T.  1)  GO  TO  100 
DO  1  1=1, N 

DO  6  J=1.N 

PI  ( I ,  J ) =CMPLX <0*0, 0,0) 

UNI T ( I , J ) =CMPL X ( 0 . 0 , 0 , 0 ) 

6  CONTINUE 

1  UNIT(I,I)=CMPLX(1, 0,0.0) 

100  CONTINUE 

****************************************************************** 

*  NONLINEAR  MULTISTEP  STARTS  HERE.  * 

*  BEGINNING  SECTION  DOES  INITIALIZATION  * 

C  ************************************* **************»******x******* 

DO  132  1  =  1 ,  N 
132  YN(I)=CMPLX(C ,0,0.0) 

IF ( IS  *  GT » 1 )  GO  TO  131 
DO  103  1=1, N 
DO  103  J=1,N 

103  P 1 ( I , J ) =-EAH ( T , J  >  +UN I T ( I , J ) 

C  ***************************************************************** « 

C  *  1ST  ORDER  GAB  * 

C  *  DO  LOOP  105  CAIC, HATES  PHI(l.O)  * 

C  *  LOOP  108  OR  11.0  COMPUTES  FINAL  Y(N+1)  * 

T  ****************************************************************** 

DO  106  1*1, N 
DO  106  J=1,N 
PH--CMPI  X(0» 0,0.0> 

DO  1 05  K>  i  , N 

105  PH=PH-A1  ( I  *K )  *P1  *K*  ' 

FE (  T  ,  l',=PH 
.1  06  CONTINUE 

*i^1  CAI  I  GFN(G,H,N,Y,l»),A,ci6*ir'* 

Tin  1  ‘l.ss’i 

DO  *•  Of* 

10A  YNfT'sYNa'+€AHI  I*  J'*Y<1«  I  )+H«FF  Cl  .  i»*H< .)) 

RETURN 
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DO  9  1*1 ,N 

DO  10  J-1,N 
C<I,J)»AA<I,J>*H/2.0 
PP<I, J)-O.DO 
10  CONTINUE 

C<I,I)«C<I, I)+1.D0 
9  CONTINUE 
DO  12  1*1 f  N 
DO  13  J*1,N 
DO  14  K*1 , N 

PP<I, J)*PP<I,J>+B<I,K>*C(K,J> 

14  CONTINUE 
13  CONTINUE 
12  CONTINUE 

IF<M  .EO.  0)  GO  TO  40 

****************************************************************** 
*  NORM < AH )  » GT « ( » 1 ) ,  EXP(A)  *EXP<A/2**M)**(2**M)  * 

****************************************************************** 

DO  24  1*1, N 

DO  25  J*1 »N 
B< I » J)»O.DO 
25  CONTINUE 

24  CONTINUE 
DO  36  K*1 , M 
DO  27  1*1, N 
DO  28  J*1,N 
DO  29  L*1,N 

B ( I , J ) *B  < I , J ) +PP  < I , L ) *PP  <  L , J  ) 

29  CONTINUE 
28  CONTINUE 
27  CONTINUE 
DO  31  1*1, N 
DO  32  J»1 , N 
BB*B< I , J) 

P<I, J)*CMPLX<BB,0.0) 

B ( I , J  >  *0 ♦ DO 
32  CONTINUE 
31  CONTINUE 
36  CONTINUE 
H-HAUE 
RETURN 
20  H»H/2 » 0 
M-M+l 

DO  54  1*1, N 
DO  55  J*1,N 
PP< I, J)*O.DO 
55  CONTINUE 
54  CONTINUE 
GO  TO  30 
40  H-HAVE 
RETURN 
END 
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START 

SUBROUTINE  START <N»T»YZERO> 

COMMON  PEPS  f  PC T  f PCN  t PBETA f  FQ  f  PA  f  PB  f  PC  f  PD » DZ  » I BND 
COMPLEX  YZERO <1> 

DIMENSION  T< 1 ) 

******  USER  SHALL  REVISE  THIS  PORTION  TO  INCORPORATE  HIS  INITIAL  VALUE 

.  YZERO<I)  CONTAINS  THE  FUNCTION 

. .  YZERO< 1+1 )  CONTAINS  THE  DERIVATIVE 

DO  29  LP“1 f Nf 2 
DLZ«DZ*(LP+l)/2 

IF  <  LP ♦ EO ♦ 1 >  YZERO < LP ) -CEXP <  CMPLX < 0 . 0 » FQ*T < 1 ) ) ) 

*  #CMPLX< . 99546739 f-.73021506E-01 ) 

IF<LP.EQ.3>  YZERO ( LP ) *CEXP < CMPLX < 0 ♦ 0 » FQ*T < 1 ) ) > 

*  *CMPLX< .98602435 14547502) 

IF(LP»E0.5)  YZERO <LP ) *CEXP<  CMPLX (0  *0  f FQ*T ( 1 ) ) ) 

*  *CMPLX< .971 72089 21705468) 

IF(LP»E0,7)  YZERO<LP)*CEXP<CMPLX<0#0f FQ*T < 1 ) ) ) 

*  *CMPLX < ♦ 96262427 » - . 28745766  > 

YZERO  <  LP+1 YZERO  <  LP ) *CMPLX  <©«0f FQ> 

29  CONTINUE 
RETURN 
END 


GFN 


SUBROUTINE  GFN <G»H»N»Y*J»T>A» SAVE ) 

C  *****  THIS  HANDLES  THE  G  VECTOR 
C  *****  G  VECTOR  CONTAINS  BOUNDARY  INFORMATION 

COMMON  PEPS  f  PCT  *  PCN f PBETA fFQf PA? PBf PC f PD fDZf  I  BND 
COMPLEX  A< 18f 18) f Y(2» 18) * G < 18 > f SURFf BOTT , XA f XBf SAVE 
DIMENSION  T ( 1 ) 

DATA  PI » ZERO > ONE/3*  1415926535f 0*0f  1  »0/ 

DO  1  1*1  fN 

1  G<I)*CMPLX<O.OfO.O> 

C  *****  FIRST  ARGUMENT  Of  CALLS  FOR  FIXED  END  CONDITION 
CALL  BC<OfNfHfTfYfSAVEfSURF> 

G  (  2  > -PC4SURF/DZ+2 . *PBETA*CMPLX ( 0 . 0  f  FQ ) *SURF/DZ 
C  *****  FIRST  ARGUMENT  If  CALLS  FOR  FREE-END  CONDITION 
CALL  BCdfNfH.TfYfSAVEfBOTT) 

X— PBETA* ( 1  ♦  - .  5*PCT*PEPS*  <  1  ♦  - < N/2  >  *DZ > )  /  <  DZ*DZ ) 

G<N) ®CMPLX(Xf 0  » 0 )*BOTT 

RETURN 

END 
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INVERT 


SUBROUTINE  INVERT<A»N»ANS) 

:M*************t************************************************** 


*  MATRIX  INVERSION  SUBROUTINE f  CALLED  BY  PADE  OR  NLMS  * 

*  A  CONTAINS  THE  ORIGINAL  ELEMENTS  AND  REMAINS  UNALTERED  * 

*  ANS  CONTAINS  THE  A**<-1)  * 

*  THIS  SET-UP  IS  USING  UNIVAC  1108  MATHPK  EXISTING  DOUBLE  * 

*  PRECISION  GAUSS- JORDAN  REDUCTION  * 

*  THIS  PROGRAM  IS  REPLACEABLE  BY  THE  USER  * 


****************************************************************** 
DOUBLE  PRECISION  A< 18» 18) t ANS<18» 18) » V<2) 

DIMENSION  JC( 18) 

DATA  NR/18/»NC/18/ 

V  < 1 ) *1 » DO 
DO  1  1-1 »N 
DO  1  J«1fN 

1  ANS  <I»J)*A(I»J) 

CALL  DG JR ( ANS  t NR  »  NC  r  N  »  N » MDEX  t JC »  V  ) 

IF<MDEX  .EG,  1)  GO  TO  10 
RETURN 

10  WRITE  (4*2) 

2  FORMAT <3X»22HMATRIX  INVERSION  ERROR) 

RETURN 

END 


INVER 


SUBROUTINE  INVER* AfNf ANS) 

PARAMETER  NDIM-18 

COMPLEX  A(NDIMtNDIM) f ANS (NDIMrNDIM) fV<2) 
DIMENSION  JC(NDIM) 

DATA  NR/NDIM/fNC/NDIM/ 
V<1)=CMPLX(1.0f0,0) 

DO  1  1  =  1  fN 
DO  1  J=1fN 

1  ANS<IfJ)=A<If J) 

CALL  CGJR(ANSfNCfNRfNfNfMDEXfJCfV) 

IF < MDEX  .EG.  1)  GO  TO  10 
RETURN 

10  WRITE  <4f2) 

2  FORMAT <3X»!.1HMAT  INV  ERR) 

RETURN 

END 
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DGJR 


SUBROUT I NE  DO JR ( A  t NC  f  NR  t N » MC » MDEX  t JC » V ) 
DIMENSION  JC<N)fV<2> 

DIMENSION  A<NR»NC) 

DOUBLE  PRECISION  ArVrXrDLOG 
IW-Wl) 

M«1 
S-l « 

L»N+<MC-N)#<IW/4> 

KD»2-  MOD( IW/2»2) 

IF<KD»EQ* 1 )  U<2>=0. 

KI=2-  MOD< IW»2) 

60  TO  (5t20) tKl 
5  DO  10  1*1 f N 

10  JCdJ-I 

20  DO  91  1*1 tH 

GO  TO  <22»21 ) »KI 

21  M-I 

22  IF  ( I . EQ.N)  GO  TO  60 

X— !♦ 

DO  30  J*I »N 

IF  < X ♦ GT * ABS <A<J»I)))  GO  TO  30 
X»ABS<A<J,I>) 

K*J 

30  CONTINUE 

IF(K.EO.I)  GO  TO  60 

mUa, 

GO  TO  <35»40) *KI 
35  MU*JC( I ) 

JC<I)-JC(K) 

JC  <  K ) *MU 

40  PO  50  J*M>L 

X*A<I, J) 

A<I»J)*A<K»J> 

50  A<K* J)»X 

60  IF  < ABS (A(IfI))«GT«Ot)  GO  TO  >0 

IF < KD  *EG . 1 )  V<1>*0. 

JC<1)*I-1 

RETURN 

70  GO  TO  <71 *72) »KD 

71  IF< A< I » I ) ♦ LT. 0» )  S*-S 

V<2)*Y<2>  +  DLOG ( ABS  <A(I»I>>> 

72  X*A< I » I ) 

A< I f I )*1  ♦ 

DO  80  J*M»L 
A<I*J)»A(I*J)/X 
CALL  ERRTST<72,MPEX> 

v  XF<MDEX.EQ.1>  GO  TO  150 

80  CONTINUE 
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DO  91  K*1»N 

IF  (K.EG.I)  GO  TO  91 
X=A<K  » I ) 

A<K,I>=0. 

DO  90  J  =M  t L 
A  <  K » J  >  =A ( K » J ) -X#A  <  I » J  ) 
CALL  ERRTST  <  72  »  MDEX ) 

IF ( MDEX . EQ ♦ 1 )  GO  TO  150 

90  CONTINUE 

91  CONTINUE 

GO  TO  <95» 140) »KI 
95  DO  130  J=1»N 

IF< JC< J) »EG» J)  GO  TO  130 
JJ*J+1 

DO  100  I*JJ»N 
IF ( JC ( I ) ♦ EQ » J)  GO  TO  110 
100  CONTINUE 
110  JC<I)*JC<J) 

DO  120  K*1 f N 
X«A<Kf I) 

A<K»I)=A<KfJ) 

120  A(K» J)=X 
130  CONTINUE 
140  JC  < 1 ) “N 

IF(KD.EQ.l)  M<1)*S 
RETURN 

150  JC<1)»1-I 

IF <KD«EQ« 1 )  V<1)*S 

RETURN 

END 
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CGJR 


SUBROUTINE  CGJR< AfNCfNRfNfMCf IFLf JCfV) 
DIMENSION  JC(1) 

COMPLEX  CL0GfUfXCfA<18f18> 

COMPLEX  Z 

INTEGER*2  NERR 

NERR-72 

IFL-0 

IU*V 

M«<0.fO.) 

IBIT-0 

M-l 

L*N+  <  MC-N  >  #  < I W/4 ) 

KD«2-M0D<IU/2f2> 

KI«2-M0D(IWf2) 

GO  TO  (5f20) fKI 
5  DO  10  I=1fN 

10  JC<I)»I 

20  DO  91  1-1, N 

GO  TO  <22f21)fKI 

21  M=I 

22  IF  < I . EQ.N)  GO  TO  60 
X— 1. 

DO  30  J»IfN 

ANORM* ABS <  REAL  <  A  <  J f 1 ) ) )  + ABS  <  A I MAG ( A( Jf I) ) ) 
IF < X « GT • ANORM >  GO  TO  30 
X«ANORM 
K= J 

30  CONTINUE 

IF<K«EQ«I)  GO  TO  60 
IBIT*IBIT+1 
GO  TO  <35f40) fKI 
35  MU»JC<I) 

JC<I)*JC<K> 

JC(K)*MU 

40  DO  50  J*MfL 

XC*A< I  f  J) 

A<IfJ)*A<Kf J) 

50  A(K»J)*XC 

60  ANORM*ABS<REAL< A< I  f  I ) ) )+ABS<AIMAG(A< I  f  I )  >  > 

IF < ANORM *GT *0)  GO  TO  70 
V»<0. fO* ) 

JC<1)»I-1 

RETURN 
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70  GO  TO  <71 f72>  fKD 

71  0=0+CL0G  <A<I»I>) 

Z*CLOG<  A< I  f  I )  > 

72  XC-A<IfI> 

A<IfI)«<1.  fO. ) 

DO  80  J«MfL 

A  ( I f  J ) “A  <  I  f  J ) /XC 
CALL  ERRTST<NERRfIFL> 
IF(IFL.EQ.l)  GO  TO  150 
80  CONTINUE 

DO  91  K»1 fN 

IF  (K.EQ.I)  GO  TO  91 
XCssA<K»  I ) 

A<KfI)»<0.f0f> 

DO  90  J  *MfL 
A(Kf J)*A(K» J)-XC*A< If J) 
CALL  ERRTST  <  NERR  f  IFL  > 
IF<IFLfEQ.1>  GO  TO  150 

90  CONTINUE 

91  CONTINUE 

GO  TO  <  95f 140 >  fKI 
95  DO  130  J=1 fN 


IF 

<JC<J)fEOfJ) 

GO 

TO 

130 

JJ= 

=J+1 

DO 

100  I=JJ»N 

IF 

< JC< I) »EG« J) 

GO 

TO 

110 

100  CONTINUE 
110  JC<I)=JC<J> 

DO  120  K=1fN 
XC=A(KfI> 

A  <  K » I ) “A  <  K  f  J ) 

120  A<KfJ)=XC 
130  CONTINUE 
140  JC<  1 ) -  N 

0=0+ <0. f 3. 14159265 >*CMFLX <FLOAT<MOD< IBITf 2) ) fO. ) 
RETURN 

150  JC<1)=1-I 

RETURN 
END 
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SUBROUTINE  BC(JrN»H»T»Y» SAVE ,  X ) 

C  *****  THIS  SUBROUTINE  SUPPLIES  THE  FREE-END  BOUNDARY  CONDITION 
C  *****  USER  HAS  AN  OPTION  TO  SUPPLY  HIS  BOUNDARY  CONDITION 
C  *****  IN  THIS  CASE?  STATEMENT  AFTER  2  SHOULD  BE  REPLACED  BY  USER'S 
C  CONDITION 

COMMON  PEPS  ,  PCT  »  PCN  *  PBETA  ,  FCJ , PA , PB ,  PC  ,  PD » D2  ,  I BND 
DIMENSION  T< 1 ) 

COMPLEX  Y(2»18)»X 
IF< J  .GT ♦  0)  60  TO  10 
C  *****  FIXED  END  BOUNDARY  CONDITION 
X=CEXP  <  CMPLX  <  0 ♦ 0  *  FQ*T ( 1 > )  ) 

RETURN 
10  CONTINUE 

GO  TO  (1,2,3) t I BND 

C  *****  R*M»KENNEDY  BOUNDARY  CONDITION 
C  ********************* 

C  *  U  +U  =0  * 

C  *  T  Z  * 

C  ********************* 

1  CONTINUE 

X- <  H*Y ( 1  *  N-l ) /DZ+SAUE )  /  <  1 , 0+H/DZ ) 

RETURN 

C  *****  USER  SUPPLIED  BOUNDARY  CONDITION 

2  CONTINUE 

X=CEXP  <  CMPLX  < 0.0,  FQ*T < 1 ) ) ) *CMPLX  < . 92986 1 96  > - . 35670799  > 
X*CMPLX(1 .0,0.0) 

RETURN 

3  CONTINUE 

IF(PBETA  .NE.  1.0)  GO  TO  31 
X=Y  < 1 »N-1 ) 

RETURN 

31  XYZ=1.+(PCT*PEPS*PBETA*DE)/(2.#<1, -PBETA) ) 

Xs® (CMPLX <  1 , /XYZ, 0 ♦  0 )  ) * < CMPLX <1.+XYZ,0.0)#Y(1, N-l  )— Y<  1  ,N-3 ) ) 

RETURN 

END 
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DIFEQ 


SUBROUTINE  DIFEQ <H»N» TMAX » Y * YZERO * ANORM t YN » SAVE ) 

COMMON  PEPS » PCT  » PCN » PBETA tFQtPAtPBtPCfPD,DZ?l BND 
COMMON  A<18f 18)fT<3) 

COMPLEX  Y<2» 18) »  YN< 18) ? Y0LD(2» 18) » YZERO < 18) » A»SAVE 
DO  40  1  =  1  »N 
40  Y  <  1  *  I  =YZERO< I ) 

IH»0 

TZERO=T  < 1 ) 

60  TEA=T ( 1 ) +H 

IF < TEA. GT. TMAX)  H=TMAX-T<1) 

IF < TEA. GT. TMAX)  IH=0 
IF < TEA. GT, TMAX)  GO  TO  60 
IH-IH+1 

IF ( IH  .GE.  32767)  IH=2 
T  <2)=T  < 1 )+H 
IMP =2 

DO  62  J=1 *  IMP 
DO  62  1*1  »N 
62  YOLD ( J  *  I ) *Y  <  J » I ) 

CALL  NLMS < H » YOLD  » NfYNrlHr SAVE ) 

59  DO  66  1*1* N 
Y<2» I )=YN< I ) 

66  Y0LD(2* I )=YN< I ) 

*jM«************************«************************************* 

*  RESULTS  Y(TEA)  IN  YN<I)  AND  Y(2*I)  * 

******************************3Mc*****************************«**** 

ANORM=TEA 

DO  85  1=1 ,N 

Y< 1 » I ) *Y < 2 *  I > 

YZERO  < I ) =Y ( 1  *  I ) 

85  CONTINUE 
T ( 1 ) =T ( 2 ) 

TZERO=T ( 1 ) 

IF < ABS ( TEA-TMAX ) ,LE . ( . IE-5) )  RETURN 

GO  TO  60 

FND 
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PADE 

SUBROUTINE  PADE< A»H»P»N> 

C  *****  A  RATIONAL  APPROXIMATION  OF  MATRIX  EXPONENTIALS 
C  *****  DOUBLE  PRECISION  IS  NEEDED  FOR  REQUIRED  ACCURACY 
PARAMETER  NM*18 
COMPLEX  A  <  NM , NM ) » P  <  NM  »  NM ) 

DOUBLE  PRECISION  AA  <  NM , NM ) » PP  <  NM » NM ) »  B  <  NM , NM  > » C  <  NM » NM ) » HAVE 
DOUBLE  PRECISION  CC  <  NM )  ,  COL  , XNORM 
HAVE=H 
DO  2  1=1, N 

DO  1  J=1,N 
B<I, J)=O.DO 
C<I, J)=O.DO 
PP( I, J)=0«D0 

AA  <  I  ,  J  ) =DBLE  <  REAL <  A  <  I  »  J  ) )  ) 

1  CONTINUE 

2  CONTINUE 
DO  17  1=1, N 
COL»O.DO 
DO  16  J»1,N 

COL=DMAXl <  COL , DABS ( DBLE  <  REAL  <  A  <  I »  J  >  >  >  >  > 

16  CONTINUE 
CC< I )=COL 

17  CONTINUE 
XNORM-CC(l) 

DO  18  1=1, N 

IF < XNORM  • 6T ♦  CC(I)>  GO  TO  18 
XNORM=CC ( I ) 

18  CONTINUE 

****************************************************************** 

*  COLUMN  NORM  IS  USED  TO  SEE  WHETHER  EXP(A)  NEEDS  REDUCTION  * 

M=0 

30  IF ( XNORM*H  -  0.98)  3,20,20 

****************************************************************** 

*  EXP (A)=(I-»  5*A )**<-l)*<I+.5*A>  * 

****************************************************************** 

3  DO  6  1=1, N 
DO  5  J=1,N 

DO  4  K=1 , N 

PP<I, J)=PP<I, J)+AA<I,K)*AA(K, J> 

4  CONTINUE 
C( I , J)«-AA< I , J)*H/2.0 

5  CONTINUE 
C<I,I)*C(I,I)+1.D0 

6  CONTINUE 
CALL  INVERT (C»N,B> 
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